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NOTATION 


A  constant  which  characterizes  the  adiabatic  nature  of  the 
liquid  medium  (for  water  B  -  3000  atmospheres) 

Isentropic  sound  speed  in  the  medium  as  a  function  of  ambient 
and  transient  pressures 

Sound  speed  in  the  undisturbed  liquid  medium 

Specific  enthalpy  of  the  liquid  medium 

A  constant  which  characterizes  the  adiabatic  nature  of  the 
liquid  medium  (for  water,  n  -  7) 

Pressure  in  the  liquid  at  the  bubble  wall 

Pressure  of  the  gas  inside  the  sphere 

Initial  pressure  of  the  gas  inside  the  sphere 

Pressure  in  the  liquid  outside  the  bubble  wall 

Pressure  in  the  undisturbed  liquid  medium;  ambient  pressure 

Instantaneous  radius  of  the  imploding  sphere 

Initial  radius  of  the  imploding  sphere 

Standoff  (measured  from  the  bubble  center);  component  in  the 
direction  of  the  radial  spherical  coordinate 

Time 

Time  measured  at  the  bubble  wall;  time  measured  when  the 
Eulerian  position  vector  r  =»  R 

Instantaneous  velocity  of  the  bubble  wall 

Eulerian  velocity  in  the  fluid  outside  the  bubble  wail 

Instantaneous  specific  volume  of  the  gas  inside  the  bubble 

Initial  specific  volume  of  the  gas  inside  the  bubble 

Polytropic  gas  constant  for  an  adiabatic  process 

Density  of  the  undisturbed  liquid  medium 
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ABSTRACT 

This  paper  presents  a  method  for  calculating  the  instantaneous  pressure, 
velocity,  acceleration,  and  radius  associated  with  the  collapse  of  a  spherical 
gas-filled  cavity  in  an  infinite  compressible  liquid.  The  method  is  an  independ¬ 
ent  approach  which  makes  use  of  Hamming's  technique  to  numerically  integrate 
Gilmore's  differential  equations  which  describe  the  collapse. 

Included  is  a  computer  program  which  will  perform  the  necessary  calcula¬ 
tions  on  a  IBM  7090/1401  digital  computer.  Results  obtained  are  in  good  agree¬ 
ment  with  those  of  Hickling  and  Plesset,  whose  work  was  unknown  to  the  present 
author  when  he  undertook  the  study. 

It  may  be  inferred  that  the  peak  shock  wave  pressure  is  significantly  re¬ 
duced  by  a  decrease  in  ambient  pressure,  an  increase  in  internal  pressure,  and/or 
a  variation  of  the  specific  heat  ratio  by  proper  selection  of  the  gas.  Control  of 
the  last  two  parameters  can  be  investigated  as  a  possible  means  of  protecting 
glass  spheres  against  sympathetic  implosion  in  multiple  sphere  buoyancy 
systems. 

ADMINISTRATIVE  INFORMATION 

This  work  was  funded  under  Special  Projects  Office  Project  Order  Number  6-0002. 

INTRODUCTION 


PURPOSE 

Because  of  the  excessive  weight-displacement  ratios  obtained  with  tough  metals  such 
as  steel  or  aluminum,  designers  are  turning  toward  nonductile  materials  for  use  in  buoyancy 
systems  for  all  depth  vehicles.  Spherical  glass  shells  are  among  the  components  for  such 
systems.1’  2  In  a  system  which  contains  a  number  of  buoyancy  spheres,  it  is  essential  to 
know  the  effect  that  the  collapse  of  one  sphere  will  have  on  neighboring  spheres  in  order  to 
prevent  catastrophic  failure.  A  two-part  investigation  has  been  initiated: 

1.  The  definition  of  the  free-fteld  pressure-time  history  due  to  the  implosion  of  a  single 
sphere. 

2.  Determination  of  the  loading  and  response  of  a  sphere  to  the  pressure  field  generated 
by  the  implosion  of  a  neighboring  sphere. 
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This  report  deals  with  the  analytical  determination  of  Part  1,  based  on  the  assump- 
tions  that  the  spherical  shell  has  negligible  weight  and  thickness  and  that  it  contains  air  at 
arbitrary  pressure. 

BACKGROUND 

The  need  for  more  complete  understanding  of  the  hydromechanical  problem  of  cavita¬ 
tion  and  the  gas  bubble  phenomena  of  underwater  explosions  has  encouraged  more  and  more 
detailed  investigations  into  the  pulsations  of  underwater  gas  bubbles.  One  of  the  earliest  of 
these  investigations  was  made  by  Rayleigh  in  1917. 3  A  more  refined  treatment,  was  success¬ 
fully  completed  by  Herring  in  1941. 4  Sometime  later  (1952),  Gilmore5  took  a  different  ap¬ 
proach  and  postulated  equations  to  describe  the  growth  or  collapse  of  a  spherical  bubble  in 
a  viscous  compressible  liquid.  Gilmore's  description  is  presented  in  the  next  section. 


THEORY 

GILMORE'S  BUBBLE  WALL  EQUATION 

On  the  basis  of  the  Kirkwood-Bethe  hypothesis,6  Gilmore  has  derived  an  equation 
(w’  ich  he  calls  a  “second  order"  approximation)  which  accurately  describes  the  (nonmigra- 
tory)  oscillations  of  a  spherical  gas-filled  cavity  in  an  infinite  compressible  liquid.  If  R  is 
the  radius  of  the  sphere,  H  the  specific  enthalpy  of  the  surrounding  liquid,  and  C  the  isen- 
tropic  sound  speed  in  the  liquid,  then  Gilmore’s  equation  is: 


M-4)-M'4M'4)4('4) 


S(0)  «  RQ,  R( 0)  -  0 


where 
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Here  cw,  end  are  respectively  sound  speed,  pressure,  and  density  in  the  undisturbed 
liquid.  B  and  n  are  constants  which  characterize  the  adiabatic  compression  of  the  liquid 
(for  water,  B  »  3000  atm,  n  ■  7).  P  is  the  pressure  in  the  liquid  at  the  bubble  wall.  If  vis* 
cosity  and  surface  tension  are  neglected  and  the  pressure  p  inside  the  bubble  is  uniform, 
then  pressure  is  continuous  across  the  boundary  of  the  sphere,  i.e.,  P  *  p  except  at  time 
t  -  0  when  the  pressure  in  the  fluid  is  artificially  and  discontinuousiy  reduced  from  P  ■  px 
to  P  •  p( 0)  «  pQ.  The  gas  can  be  assumed  to  undergo  an  adiabatic  expansion  (or  compression). 
From  thermodynamics,  for  an  ideal  gas, 


P0»X  “  P*y 

where  y  is  a  constant  (the  speciric  heat  ratio),  «  is  specific  volume,  and  the  subscript  0  refers 
to  some  initial  state.  Since  the  volume  of  a  sphere  is  proportional  to  its  radius 
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Elimination  of  volumq  between  Equation  [1.4a]  and  Equation  [1.4b]  yields 


P  -  < 


Pm 


*n\  3y 


\pmp°  (t) 


t  m  0 


t  >  0 


[1.4c] 


For  air  (y  -  4/3),  the  exponent  3y  becomes  4.  Usually  the  value  of  y  is  taken  to  be  1.4  for 

air.  This  value  represents  the  behavior  of  air  fairly  accurately,  but  since  y  decreases  with 
•  •  ' 

increasing  pressure,  4/8  represents  a  rough  average.  It  will  be  seen  later  that  the  use  of  a 
constant  value  of  y  leads  to  a  deficiency  in  the  model. 
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Combining  Equations  [1.1],  [1.2],  [1.3],  and  [1.4c]  and  taking  n  -  7  yields  the  folio- 
ordinary  differential  equation  for  R : 


R(0)  -  R0,  R( 0)  -  0 


Once  B,  RQ,  pQ ,  pm,  and  cm  we  specified,  it  is  possible  to  find  a  numerical  solution  for 
R(t),  R(t),  R(t),  and  p(t). 

The  initial  velocity  i?( 0)  is  taken  throughout  this  paper  to  be  zero.  With  the  help  of 
Equation  [1.5],  Oilmore  has  pointed  out  that  near  t  ■  0,  there  is  a  small  finite  jump  in  velocity 
during  an  infinitesimally  small  interval  of  time,  i.e.. 


*(0+) 


Pq-Po. 

PmCm 


Hickling  and  Plesset7  give  a  good  physical  explanation  of  this  jump  in  terms  of  the  initial 
pressure  discontinuity  between  pQ  and  This  velocity  jump  may  lead  cae  to  choose 
*<0+)  -  (pQ  -  p„)/patctm  as  the  initial  condition  on  the  velocity.  Since  the  difference 
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between  (pQ  -  Pm)/pm°m  And  zero  is  small  compared  to  the  magnitude  of  the  velocities  of 
interest,  the  question  as  to  whether  to  start  the  solution  at  i?(0)  ■  0  or  at  R(Q^)  »  (p0  -  pm)/ 
p  em  is  somewhat  academic.  Tho  plots  of  bubble  wall  velocity  do  not  show  this  initial  jump 
because  the  writer's  solution  of  Gilmore's  equation  was  actually  carried  out  from  t  -  0+  to 
avoid  an  infinite  initial  acceleration.  (One  of  the  terms  appearing  in  the  expression  for  the 
initial  acceleration  is  the  derivative  with  respect  to  time  of  P  as  given  in  Equation  [1.4c]; 
this  derivative  is  infinite  at  t  •  0).  The  approximation  (p0  -  pj)/pmcm  m  0  was  made  to 
simplify  the  calculations  of  the  initial  values  of  R  and  R  which  are  tedious  even  with  such 
an  approximation. 

It  shbuld  be  emphasized  that  one  of  the  inherent  assumptions  upon  which  Equation 
[1.5]  is  based  is  that  the  cavity  remain  spherical  throughout  the  collapse  and  subsequent 
oscillations.  In  some  bubble  collapse  experiments,8  however,  the  single  bubble  has  occa¬ 
sionally  been  observed  to  dissociate  into  many  smaller  bubbles  at  the  end  of  the  first  collapse. 

THE  INTEGRATION:  HAMMING’S  METHOD 

An  ordinary  differential  equation  of  the  form 

—  y(*0)«y  o  [2.1] 

can  be  integrated  numerically  by  one  of  various  finite  difference  methods.  One  of  these,  the 
Hamming  method,  which  is  particularly  well  suited  for  the  solution  of  Equation  [1.5]  can  be 
found  in  Ralston  and  Wilf.9  It  is  outlined  briefly  here: 

1.  The  *-axis  is  equally  divided  into  a  large  number  of  small  intervals.  The  value  of  y 
at  the  end  of  the  nth  interval  (i.e.,  the  interval  between  xn  and  sn-1)  is  denoted  by  y#. 

2.  Knowing  the  values  of  y.  and  y ^  at  previous  intervals  xt  up  to  and  including  the  nth 
(f  -  n),  it  is  possible  to  calculate  p„+1,  a  first  approximation  to  the  (n  +  l)#t  value  of  yn 
at  *B+1,  by  means  of 


-  y„_3  ♦  j  (2 y;  -  +  2 y;_2) 

whose  h  is  the  width  of  the  interval;  y(l+j  is  called  the  predictor. 

8.  Since  the  prediction  pB+1  is  based  on  the  value  of  a  series  expansion  of  y,  error  due 
to  truncation  is  incurred.  Most  of  the  difference  between  the  true  value  of  y  and  the  estimated 
value  of  y  is  taken,  into  account  by  the  modifier  (denoted  by  «n„+1) 
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Wi»+  l*?n+l“  j21  ^n~Cn) 

where  on  is  as  defined  in  the  next  paragraph  and  the  derivative  of  is  m'+1  * 
f(Xn+l’  mn+ 1)* 

4.  The  predicted  value  is  compared  with  a  quantity  called  the  corrector 

c»+i  *  \  1 +  3A  «+i +  2y»  - 

5.  If  the  predictor  pn+1  lies  close  to  cn+i  within  some  specified  tolerance,  then  the  final 
value  of  yn+1  at  xn+l  is  taken  as 


^n+l  ”  Cn+1  +  ^21  ^n+t  Cb+1^ 

6.  If  the  predictor  pn+l  does  not  lie  close  enough  to  cn+1>  then  either  (1)  a  new  value  of 


9 

^n+1  "  Cb+1  +  j^21  ^ ~  ^a+l' 

may  be  calculated  and  an  interative  process  carried  out  or  (2)  the  interval  may  be  halved. 
This  is  discussed  in  more  detail  in  Appendix  A. 

In  order  to  solve  Equation  [1.5]  by  the  procedure  just,  outlined,  it  is  necessary  to 
reduce  the  equation  to  a  system  of  two  simultaneous  differential  equations  ’of  the  form  of 
Equation  [2.1].  This  method  is  mentioned  in  Hildebrand. 10  Write  V  for  R\  then  Equation 
[1.5]  can  be  written 


k  -  U 


[1.6a] 


A 


and 


Hamming’s  method  can  be  applied  simultaneously  to  Equations  [1.6a]  and  f  1.6b]  to  yield  a 
numerical  solution.  The  justification  for  use  of  this  particular  method  is  du  cussed  in 
Appendix  B. 

The  function  R(t)  (from  which  U(t)  and  p(t)  at  the  bubble  wall  can  be  obtained),  deter¬ 
mined  by  Equation  [1.6],  constitutes  one  of  the  boundary  conditions  necessary  to  find  the 
Eulerian  velocity  and  pressure  fields  in  the  fluid  outside  the  bubble  wall. 


THE  EULERIAN  VELOCITY  AND  PRESSURE  FIELDS  IN  THE  LIQUID 

In  his  “second  order"  approximation,  Gilmore  uses  the  Kirkwood-Bethe  hypothesis  in 
conjunction  with  the  method  of  characteristics  to  determine  the  (Eulerian)  velocity  and  pres¬ 
sure  fields  in  the  liquid.  If  the  standoff  r  is  greater  than  or  equal  to  the  initial  radius  RQ, 
then  the  velocity  u  associated  with  the  standoff  will  not  be  of  the  same  order  of  magnitude 
as  the  sound  speed  except  for  the  most  severe  implosions.  Provided  the  approximation 
t»2  «  c2  is  valid,  the  following  set  of  equations  (the  expressions  derived  by  Gilmore)  are 
sufficient  to  determine  the  velocity  and  pressure  fields  u  and  p  in  the  liquid  when  U  and  R, 
the  bobble  wall  velocity  and  radius,  are  known  functions  of  time. 


CJf  f  / 


[3.1] 


7 


where  j r  and  K3  are  given  by 


Here  tR  is  the  time  at  which  the  bubble  radius  is  R  and  the  bubble  wall  velocity  is  U.  An 
event  which  occurs  at  the  bubble  wall  at  time  tR  requires  finite  time  t  -  tR  to  propagate  Irom 
the  bubble  wall  to  the  point  r  in  the  fluid.  This  time  lag  is  indicated  by  Equation  [3.5]. 

If  the  approximation  u2  «  c 2  is  not  made,  then  the  numerical  integration  for  u  and  p 
is  more  complicated  and  requires  a  great  deal  more  time  on  the  computer.  Hickling  and 
Plesset?  avoid  making  this  approximation.  Instead  of  using  values  of  U  and  R  in  Equations 
[8.1]  through  [3.5]  to  find  p  and  u  at  discrete  points,  they  use  V  and  R  as  coefficients  in  a 
differential  equation* which  determines  p  and  v.  Then  each  time  U  and  R  are  determined  at 
a  single  point,  another  differential  equation  must  be  solved  to  find  p  and  u.  Each  solution 
of  this  second  differential  equation  gives  p  and  u  as  a  function  of  distance  from  the  bubble 
at  one  specific  instant  in  time  (i.e.,  that  instant  in  time  at  which  the  bubble  radius  is  that 
value  of  R(t)  used  in  this  second  differential  equation). 


RESULTS 

The  integration  of  Equation  [1.8]  (the  equation  which  determines  bubble  radius  and 
bubble  wall  velocity  and  acceleration) "has  been  coded  in  FORTRAN  for  a  7090/1401  com¬ 
puter  according  to  the  procedure  outlined.  The  program  and  some  sample  input  and  output 
are  given  in  Appendix  C.  Plots  from  computer  output  of  R,  u,  0,  and  p  (the  bubble  radius, 
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bubble  wall  velocity  and  acceleration,  and  pressure  at  the  bubble  wall,  respectively)  as 
functions  of  time  can  be  found  among  Figures  1  through  5  for  various  ambient  and  internal 
pressures. 

Equations  [3.1]  through  [3.5]  (those  equations  which  determine  the  Euierian  velocity 
and  pressure  fields  in  the  fluid  outside  the  bubble  wail)  have  been  incorporated  into  the 
program.  From  computer  output,  plots  of  u  and  p  (Euierian  velocity  an  ressure)  were  ob¬ 
tained  and  are  found  with  the  corresponding  plots  of  R ,  U,  U,  and  p  (Figures  1  through  4). 

DISCUSSION 

Dynamically,  the  air  inside  the  bubble  behaves  in  a  peculiar  fashion  which  is  quite 
evident  in  the  more  violent  implosions  (those  at  great  depths).  During  the  greater  part  of 
the  collapse,  the  air  offers  insignificant  resistance  to  the  inrushing  water.  Just  before  the 
instant  of  minimum  radius,  however,  the  air  violently  arrests  further  decrease  in  volume, 
behaving  very  much  like  a  rigid  sphere.  This  is  borne  out  especially  by  the  curves  for  bub¬ 
ble  wall  acceleration.  At  the  instant  of  minimum  radius,  the  water  in  the  immediate  vicinity 
of  the  bubble  “sees'*  a  rigid  sphere,  but  that  water  a  little  further  from  the  wall  continues 
to  rush  in  since  the  water  is  compressible.  The  result  is  a  spherical  shock  wave  propagating 
from  the  wall  out  into  the  fluid. 

Comparisons  among  Figures  1  through  5  lead  to  the  following  observations: 

1.  As  the  depth  at  which  the  implosion  occurs  becomes  greater,  the  peak  pressure  in¬ 
creases,  the  rise  time  decreases,  and  the  collapse  time  decreases. 

2.  Increasing  the  initial  internal  pressure  of  the  gas  inside  the  sphere  has  roughly  the 
same  effect  as  decreasing  the  depth  of  implosion.  Specifically,  the  peak  pressure  can  be 
effectively  attenuated  by  increasing  the  initial  internal  pressure.  Comparison  between  Fig¬ 
ures  2e  and  3e,  for  example,  shows  that  the  peak  pressure  pulse  from  an  implosion  at  1000  ft 
of  water  is  reduced  by  almost  40  percent  when  the  initial  internal  pressure  is  increased  from 
1  to  2  atm. 

3.  As  the  pressure  peak  propagates  away  from  the  cavity  wail,  it  suffers  an  attenuation 
proportional  to  1/r. 

4.  Comparison  between  Figures  1,  2,  and  5  indicates  that  for  constant  initial  internal 
pressure  and  radius,  the  collapse  time  varies  approximately  inversely  with  the  square  root 
of  the  depth  at  which  the  implosion  occurs. 

It  should  be  noted  that  the  present  theory  does  not  account  for  the  effect  of  migration. 
The  model  presented  here  is  excellent  for  a  ratio  of  initial  internal  pressure  to  ambient  pres¬ 
sure  which  is  less  than  perhaps  one-tenth. 
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Figure  '4  —  Spherical  Collapse  as  a  Function  of  Time  for  a 
Water  Depth  of  1000  Feet,  an  Initial  Radius  of  1  Inch, 
and  an  Initial  Internal  Pressure  of  10  Atmospheres 
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Figure  5  -  Spherical  Collapse  as  a  Function  of  Time  for  a 
F Water  Depth  of  10,000  Feet,  an  Initial  Radius  of  1  Inc  , 
and  an  Initial  Internal  Pressure  of  1  Atmosphere 
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Figure  Sb  -  Bubble  Wall  Velocity 
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Figure  5c  -  Bubble  Wall  Acceleration 
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Figure  5d  —  Preaaure  at  the  Bubble  Veil 
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Figure  St  -  (Eulerian)  Velocity 
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As  stated  earlier,  the  “well  behaved"  spherical  collapse  of  the  model  occasionally 
may  not  conform  to  the  behavior  of  a  real  bubble  in  its  final  stage  of  implosion.  At  very 
high  ratios  of  ambient  to  initial  internal  pressure,  the  possible  dissociation  of  the  cavity 
into  numerous  smaller  bubbles  represents  a  departure  from  the  behavior  of  the  model.  At  the 
standoff  of  interest  (r  >  RQ),  however,  the  field  variables  (pressure  and  Eulerian  velocity) 
are  thought  not  to  deviate  significantly  from  values  obtained  using  the  model. 

All  the  results  mentioned  so  far  may  be  extended  to  cases  for  spheres  of  any  radius. 
Suppose  that  at  depth  A  a  solution  exists  for  a  sphere  with  initial  radius  RQ  and  initial  inter¬ 
nal  pressure  pQ.  The  radius,  velocity,  acceleration,  and  pressure  are  known  functions  of 
time  at  the  bubble  wall  or  at  some  standoff  in  the  fluid.  If  the  initial  radius  is  multiplied  by 
A  -  constant,  then  pressure  and  velocity  will  remain  the  same  if  radius,  time,  and  standoff 
are  multiplied  by  A  and  acceleration  divided  by  A. 

After  the  writer’s  program  was  completed,  it  was  discovered  that  Hickling  and  Plesset 
had  solved  the  free-field  implosion  problem  numerically  in  a  similar  but  more  elaborate  man¬ 
ner  without  making  the  approximation  v2  «  c2.  Their  program  requires  20  min  of  computer 
time  for  each  case  compared  to  only  2  min  for  the  program  presented  here.  They  report  two 
cases;  the  first  (pQ  ■  10“3  atm,  px  ■  1  atm)  was  solved  by  the  program  presented  here,  but 
the  second  (pQ  *  10“4  atm,  =»  1  atm)  is  too  violent  an  implosion  for  the  writer’s  program 
to  be  applicable  because  the  approximation  u2  «  c2  may  not  be  valid.  (Note  that  this  sec¬ 
ond  case  reported  by  Hickling  and  Plesset  represents  such  a  violent  implosion  that  it  has 
little  in  common  with  the  type  of  implosion  expected  in  a  buoyancy  sphere  system  even  as 
deep  as  30,000  feet  of  water). 

In  an  attempt  to  verify  the  soundness  of  his  approach,  the  writer  used  the  initial  con¬ 
ditions  of  Hickling  and  Plesset’s  first  case  (pQ  •  10~3  atm,  =»  1  atm)  as  input  in  his  pro¬ 
gram.  Very  little  discrepancy  (less  than  2  percent)  can  be  found  in  the  bubble  radius  and 
bubble  wall  velocity  even  though  the  Hickling  and  Plesset  results  are  based  on  a  value  of 
y  ■  1.4  and  the  writer's  are  based  on  a  value  of  y  *  4/3.  It  can  be  seen  from  the  differential 
(Equation  [1.51)  that  a  small  change  in  y  has  little  effect  on  the  radius-time  curve. 

However,  a  marked  difference  appeared  between  the  results  of  the  two  programs  when 
the  peak  pressures  inside  and  outside  the  bubble  were  compared.  The  Hickling  and  Plesset 
results  showed  peak  pressures  which  were  about  twice  those  obtained  in  this  study.  This 
discrepancy  can  be  readily  resolved  by  noting  the  different  values  of  y  used.  The  results  of 
both  programs  indicate  that  a  minimum  radius  *  0.0170  will  be  obtained  when  a  sphere 
of  initial  radius  RQ  -  1  and  initial  internal  pressure  pQ  -  10“3  atmospheres  is  imploded  at 
the  ambient  pressure  pm  -  1  atm.  The  pressure  at  the  boundary,  by  Equation  [1.4c],  is 


When  y  -  1.4 


« 

?MAX  ”  ?0 


and  when  y  -  4/3 


l-y.x-fo^)4'0 


[4.1] 


[4.2] 


By  dividing  Equation  [4.1]  by  Equation  [4.2],  PUAX  is  (1/.0170)'2  or  .1.25  times  PHAX •  Had 
a  value  of  y  -  1.4  been  used  in  the  writer's  program,  the  peak  pressure  at  the  bubble  wall 
would  have  compared  well  with  that  obtained  by  Hickling  and  Plessel.  A  similar  statement 
is  true  for  those  peak  pressures  in  the  fluid  outside  the  bubble  wall,  because  the  peak  pres¬ 
sure  varies  inversely  as  the  distance  from  the  center  of  the  bubble  (i.e.  as  1/r).  This  veri¬ 
fies  the  validity  of  the  writer’s  program  and  the  assumption  u2  «  c2  for  the  range  of  interest 
[Po.  “  1000  atm,  pQ  -  1  atm). 

The  verification  against  the  work  of  Hickl :  %  and  Plesset  shows  that  small  variations 
in  y,  the  specific  heat  ratio,  can  lead  to  large  variations  in  the  peak  pressure  associated 
with  a  collapse.  Such  behavior  suggests  that  the  present  equation  of  state,  the  ideal  gas 
law,  is  a  deficient  description  of  the  gas  inside  the  cavity  and  that  the  use  of  a  more  elabo¬ 
rate  equation  of  state  (e.g.,  the  Beattie-Bridgeman1 1  equation  of  state)  would  give  more  ac¬ 
curate  results.  Use  of  the  Beattie-Bridgeman  equation  could  be  made  to  investigate  the  dif¬ 
ferences  in  behavior  of  the  collapse  for  different  gases. (representing  different  values  of  y). 
Variations  of  the  kind  of  gas  inside  the  cavity  along  with  variations  in  its  initial  internal 
pressure  may  be  used  to  control  the  characteristics  of  the  pressure  pulse  emitted  when  a 
glass  buoyancy  sphere  collapses.  Such  control  might  ultimately  be  used  to  reduce  the  dis¬ 
tances  between  glass  spheres  in  buoyancy  sphere  systems  without  increased  risk  of  sympa¬ 
thetic  implosions. 


SUMMARY  AND  CONCLUSIONS 

1.  A  program  to  integrate  Gilmore's  equations  describing  the  collapse  of  a  rpherical  gas 
filled  cavity  has  been  written  (see  Appendix  C).  The  program  is  general  enough  to  be  used 
in  the  study  of  such  phenomena  as  cavitation  and  underwater  explosion  gas  bubble  pulses, 
provided  behavior  is  adiabatic. 

2.  Parameters  from  the  program  for  various  ambient  and  initial  internal  pressures  have 
been  plotted  (see  Figures  1  through  5). 


8.  The  pressure  shock  wave  associated  with  the  implosion  may  be  controllable  in  two 
different  ways: 

a.  Proper  variation  of  the  initial  internal  pressure  of  the  gas  inside  the  sphere. 

b.  Proper  variation  of  the  specific  heat  ratio  y  by  changing  the  kind  of  gas  in¬ 
side  the  sphere. 

These  two  effects  should  not  be  overlooked  as  possible  means  of  protecting  glass  buoyancy 
spheres  from  sympathetic  implosions. 

FUTURE  WORK  PLANNED 

1.  An  experimental  verification  will  be  carried  out  to  determine  how  good  Gilmore’s 
model  is. 

2.  The  velocity  and  pressure  curves  can  be  used  to  synthesize  analytical  functions  to 
describe  the  free-field  implosion.  This  information  will  then  form  a  foundation  for  the  anal¬ 
ysis  of  the  effects  of  a  single  implosion  in  a  system  of  buoyancy  spheres.  Such  an  analysis 
has,  in  fact,  been  started. 

3.  The  present  computer  program  is  now  being  altered  by  replacing  the  ideal  gas  law  by 
the  Beattie-,  idgeman  equation  of  state. 
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APPENDIX  A 

A  PROCEDURE  FOR  HALVING  THE  INTERVAL  OF  INTEGRATION 


The  description  of  Hamming’s  method  pointed  out  that  if  the  predictor  does  not  lie 
close  enough  to  the  corrector,  then  one  alternative  is  to  halve  the  interval.  According  to 
Ralston  and  Wilf,®  a  suitable  set  of  interpolation  formulas  for  Hamming's  method  is: 

s'.-w  •  (“s',  *  135y,-i  *  wy.-j  *  »„-3)  *  555  * 

*  i35y._,  ♦  *  y„.,)  ♦  ± 

where  h  is  the  original  step  size. 

In  the  particular  program  which  was  written  for  the  solution  of  Equation  [1.5],  the 
following  procedure  was  adopted: 

1.  If  pn+1  was  not  close  enough  to  cn+1  then  an  iteration  was  carried  out. 

2.  If,  after  the  first  iteration,  the  value  of  pn+1  was  still  not  clcse  enough  to  c(|+1,  then 
the  interval  was  halved  and  the  entire  integration  process  at  that  step  was  carried  out  from 
the  beginning  using  the  new  half  interval. 

As  the  radius  of  the  bubble  approaches  its  minimum,  the  radius  and  wall  velocity 
become  more  and  more  difficult  to  predict,  that  is,  it  becomes  more  and  more  likely  that 
Pn+l  will  not  fall  within  the  desired  limits  of  cn+1>  With  the  above  procedure,  more  coordi¬ 
nates  will  be  calculated  near  the  minimum  where  the  functions  are  changing  most  rapidly 
than  are  calculated  in  regions  of  small  slopes. 


APPENDIX  B 

JUSTIFICATION  FOR  THE  USE  OF  HAMMING'S  METHOD 


The  choice  between  an  elaborate  procedure  like  Hamming's  method  and  some  other 
iterative  routine  to  integrate  the  equations  is  easy  to  make  if  it  is  based  on  economy  of  com* 
puter  time.  When  an  iterative  technique  is  used  to  converge  on  the  correct  value  of  the  de* 
pendent  variable  at  each  step,  the  finite  difference  form  of  the  differential  equation  may  have 
to  be  evaluated  many  timec  because  the  initial  prediction  of  the  dependent  variable  is  not 
likely  to  be  very  accurate.  In  the  case  of  Equations  [1.5],  it  is  desirable  to  avoid  the  eval* 
uation  of  the  finite  difference  equation  as  often  as  possible  because  it  involves  so  many 
calculations!  Hamming's  method,  on  the  other  hand,  makes  a  much  more  accurate  initial  pre* 
diction  of  the  dependent  variable  at  each  step  simply  because  more  information  about  past 
values  is  utilized  in  making  such  predictions.  Consequently,  Hamming’s  method  yields  a 
more  rapid  convergence  because  the  finite  difference  form  of  the  equation  has  to  be  evaluated 
only  once  or  twice  at  each  step  to  obtain  an  accurate  value  of  the  dependent  variable  there. 

One  disadvantage  in  using  this  technique  is  that  it  is  not  self-starting.  Values  of  R, 

V ,  and  N  in  at  least  four  equally  spaced  intervals  near  t  •  0  are  required.  Such  values  may 
be  computed  by  expressing  R(t)  in  a  Taylor  Series  in  t  about  zero  and  using  Equations  [1.5] 
to  evaluate  the  coefficients.  This  approach  for  solving  differential  equations  can  be  found 
in  any  elementary  text  on  the  subject,  for  instance,  Coddington.12  Because  of  the  difficulty 
of  differentiating  Equation  [1.5],  Herring’s4  equation 
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rather  than  Equation  [1.5]  was  used  to  find  the  first  four  values  of  R ,  (/,  and  V.  This  equa¬ 
tion  has  also  been  derived  by  Gilmore  and  called  the  "first  order  approximation." 

As  mentioned  ‘n  the  discussion  under  Equation  [1.5],  to  find  the  initial  acceleration  of 
the  bubble  wall,  Equation  [5.1]  should  be  evaluated  at  t  ■  0+  rather  than  at  t  »  0  in  order  to 


avoid  an  infinite  initial  acceleration.  It  can  be  seen  from  Equation  [1.4c]  that  Equation 


[5.1]  evaluated  at  t 


0  contains  the  infinite  term 


d2R  Pp-Pm 

" &  '  <,-c- 


.  The  result  is 
o 


in  agreement  with  acoustic  theory. 


APPENDIX  C 

THE  COMPUTER  PROGRAM  FOR  IBM  7090 


The  following  is  a  list  of  FORTRAN  IV  symbols  used  in  the  program. 


FORTRAN  Symbol 

Corresponding 
Symbol  Used 
in  Discussion 

Explanation 

A 

R 

Instantaneous  radius  of  the  imploding  sphere. 

B 

B 

A  constant  which  characterizes  the  adiabatic  na¬ 
ture  of  the  liquid  medium  (for  water,  B  «  3000  atm) 

BDYP 

P 

Pressure  in  the  liquid  at  the  bubble  wall 

B0 

*0 

Initial  radius  of  the  imploding  sphere 

B0D 

- 

A  dimensioned  variable  name  under  which  the  B0, 
the  unitial  radius  is  read  in 

B2,  B3,  B4,  B5,  B6 

- 

Coefficients  of  the  Taylor  Series  expansion  of  R 
about  t  ■  0 

C 

Coo 

Sound  speed  in  the  undisturbed  liquid  medium 

D 

Poo 

Density  of  the  undisturbed  liquid  medium  (for 
water,  px  »  2  slugs/cu  ft) 

DU 

• 

u 

Instantaneous  acceleration  of  the  bubble  wall 

H 

— 

Depth  at  which  the  implosion  takes  place 

HD 

— 

A  dimensioned  variable  name  under  which  H,  the 
depth  of  implosion,  is  read  in 

PA 

- 

Atmospheric  pressure  (14.7  psi) 

PL 

Poo 

Pressure  in  the  undisturbed  liquid  medium 

P0 

Po 

Initial  pressure  of  the  gas  inside  the  sphere 

P0D 

— 

A  dimensioned  variable  name  under  which  P0,  the 
initial  internal  pressure,  is  read  in 

R  and  T 

*R 

Time  measured  at  the  bubble  wall 

S 

h 

The  step  size,  or  size  of  the  interval  over  which 
the  integration  is  to  be  performed 

STF1,  STF2, 

STF(I)  (I  -  1,  2,  3) 

The  five  standoffs  for  which  Eulerian  velocities  and 
pressures  are  calculated;  note  STF(l)  =  STF3,  etc. 

STFD1, 

STFD2,  STFD(I) 

(I  -  1,  2,  3) 

Dimensioned  variable  names  under  which  the  above 
five  standoffs  are  read  in 

FORTRAN  Symbol 

Corresponding 
Symbol  Used 
in  Discussion 

Explanation 

STFPC1,  STFCP2, 
STFPC(I) 

I  -  1,  2,  8 

7 

Five  instantaneous  pressures  evaluated  at  the 
standoffs  STF1,  STF2,  STF3,  STF4  and  STF5, 
respectively;  note:  STFPC(l)  •  STFPC3  etc. 

STFUC1,  STFUC2, 
STFUC(I) 

I  -  1,  2,  8 

u 

Five  instantaneous  Eulerian  velocities  evaluated 
at  the  standoffs  STF1,  STF2,  STF3,  STF4,  and 
STF5,  respectively 

T  and  R 

*R 

Time  measured  at  the  bubble  wall 

TSTF1,  TSTF2, 
TSTF(I) 

I  -  1,  2,  8 

t 

The  time  for  each  of  the  five  pressures  and  veloc¬ 
ities  evaluated  at  each  of  the  five  standoffs, 
respectively 

U 

V 

Instantaneous  velocity  of  the  bubble  wall 

Y 

y 

A  constant  used  in  the  expression  for  the  Eulerian 
velocity  (see  Equations  [3.1]  and  [3.2]) 

YK 

Kz 

A  constant  used  in  the  expression  for  the  Eulerian 
velocity  (see  Equations  [3.1]  and  [3.3]) 

AP4,  AP5 

JV  pn+i 

Predicted  values  of  the  radius,  AP  5  is  the  pre¬ 
dicted  value  being  tested  at  the  (n  +  1)  th  interval, 
and  AP4  is  the  previous  predicted  value  of  R 
which  was  cloest  to  the  actual  value  of  R  at  the 
nth  interval 

A4MH,  A3MH 

y n—l/2'  Yn-3/2 

When  R  at  the  (n  +  1)  th  interval  is  being  predicted 
and  the  half  interval  routine  is  required,  then  these 
•re  the  interpolated  values  of  R  between  the  nth 
and  (n  -  1)  th  interval  and  between  the  (n  -  1)  th 
and  (n  -  2)  th  interval,  respectively 

C0D 

"n+l*  m»+l 

Modifier  for  U,  also  derivative  of  the  modifier  for 

R 

C4,  C5 

°e»  C„+l 

Correctors  for  l)\  C5  is  the  corrector  being  tested 
at  the  (n  +  l)th  interval  and  C4  is  the  corrector  at 
the  previous  interval 

D0D 

m»+l 

Modifier  for  R 

DC0D 

Derivative  of  the  modifier  for  U 

D4,  D5 

<V  ^+t 

Correctors  fcr  R  (see  C4,  C5) 

UP4,  UP5 

?n+l 

Predicted  values  of  the  velocity  (see  AP4,  AP5) 

U4MH,  U3MH 

^ n—l/2’  ^11-3/2 

Half  interval  values  of  velocity  (see  A4MH,  A3MH) 
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THE  COMPUTER  PROGRAM 

The  simplified  flow  chart  is  shown  in  Figure  6. 


START 


Figure  6  -  Simplified  Flow  Chart 
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DATA  INPUT 

The  first  data  card  is  read  according  to  the  format  (12)  and  must  have  a  number  no 
greater  than  50  in  the  first  two  columns  (see  Figure  7).  This  number  should  equal  the  num¬ 
ber  of  data  cards  to  follow.  Each  of  the  following  data  cards  should  contain  the  initial  infor¬ 
mation  for  a  single  collapse.  Eight  pieces  of  information  are  placed  on  each  card;  the  card 
is  divided  equa'ly  into  eight  parts,  10  spaces  each.  Information  is  read  in  according  to  the 
format  (6F10.4).  The  first  two  pieces  of  information  are  the  depth  of  collapse  in  feet  of 
water  and  the  initial  bubble  radius  in  inches,  respectively.  The  next  five  divisions  are  set 
off  for  five  values  of  standoff  in  inches  from  the  bubble  center  at  which  a  pressure  and 
(Eulerian)  velocity  time  history  are  desired.  The  last  division  is  reserved  for  the  initial 
internal  pressure  inside  the  bubble  in  pounds  per  square  inch.  The  first  eight  pieces  of 
information,  corresponding  to  one  collapse,  will  require  no  more  than  2  min  of  running  time. 
Each  additional  card,  i.e.,  each  additional  collapse,  requires  no  more  than  1  min. 

DATA  OUTPUT 

For  each  data  card,  the  computer  will  print  2000  lines  of  output.  The  first  1000  lines 
of  numbers  are  printed  under  headings  TIME,  RADIUS,  VELOCITY,  ACCELERATION,  BDYP, 
TSTF1,  STFPC1,  STFUC1,  TSTF2,  STFPC2,  STFUC2  (see  Figure  8).  The  first  five  quan¬ 
tities  refer  to  the  bubble  wall;  BDYP  is  the  internal  pressure  inside  the  bubble  (absolute 
pressure,  not  overpressure).  Each  line  refers  to  the  state  of  the  motion  at  one  instant  in 
time.  STFPC1  and  STFUC1  are  the  overpressure  (pressure  above  ambient  pressure)  and 
(Eulerian)  velocity  as  functions  of  the  time  TSTF1  for  the  first  standoff  given  in  the  data. 
Similarly  for  STFPC2,  STFUC2,  -and  TSTF2.  The  next  1000  lines  give  similar  information 
for  the  third,  fourth,  and  fifth  standoffs  specified.  The  headings  are  TSTF3,  STFPC3, 
STFUC3,  TSTF4,  STFPC4,  STFUC4,  TSTF5,  STFPC5,  STFUC5.  All  time  is  given  in  milli¬ 
seconds,  the  radius  in  inches,  all  velocities  in  inches  per  second,  the  acceleration  in  inches 
per  second  per  second,  and  all  pressures  in  pounds  per  square  inch. 
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Figure  7  -  Data  Input  for  Computer  Program 


INNNNNN«il<M)lMINItNNM«INNNIlMN«)NNI|NANnMARAMnnillMnnM)nMn*)nnMnnHI 

I00090000000000000000000000000000000000900000900000000004 


5MMMMMWWUUO0WMWW _ _ _ „ _ _ _ _  __  _ _ _ _ 

uO09®kooo  —  00000009090000900000  —  00*00  —  0990000000000900990000 

fcaoo09O0090O900o— 900000-9  9*-?99noooo-------<«  00000000000000- 

—  —  —  -.S  —  —  —  —  —  —  —  — 


IMiVUII««iMMMMMNil|ll 


m  v  w  w  i 

=  3:i: 


ivwiMMiMaiaitl 


IIHUIIVI 

0  c  o  o  0  1 

4  9  o  *.  *  < 


•  •••••lit 

m  Aunnnit  nm 

••••••••«< 

•S3SS2253S! 

tS35S5Ss::s: 

.MMMMNNMM  MM  I 


I5i?23 

1  —  M  «  • 


13383: 

>  n  •  n  •  < 

'  «  NON  1 

>*9004 


1SSS33! 

I  0  K  O  O  0  • 

»  —  —  ft  9  O  4 


IX?SS3 

•  O  O  9  0  0 


9  0  0—9 
00000 

NKNS  H 


W  II  II  M 
O  -  0  O 
#000 

9  0  0  0 

-♦to 


oeeoooooo  < 

i  i  i  r  •  •  •  i  i 

iOKOaoOKOO 
o«*«»n«A« 

►Goooconec 


»  o  o  o  0 

I  9  0  0  <0 

>  o  o  -  « 


o«  o  •  i 

-»e  ♦  i 
0  k  0  o  i 
o  o  0  9  . 
0  0  0  0  < 

NANAI 


i  M  N  I  M 
l  9  0  9  0  0  4 

io4«Ml>i 
mON«0AI 
I  «  K  K  9  9  4 

INNNNNI 


I  O  K  9  0  » 
•  *  0  O  9  A- 
f  O  9  0  K  0 

►  —  —  0  0  0 
\  0  0  0  0  0 


0>O  «  xQ 
9  0  0  0  0 
9  A«*N  A 
0-00- 


-  0  -  9 
0  0  9  9 
9  0  —  K 

K  *  —  K 
K  O  9  9 

**>  —  n  0 


<*  no  jo  o  o  i 

|5?S25238S 

ftSSSJSSc 


tOOQO  9009000001 


>0090  0  00000  00000090000000000000000  i 


llllllllllllllltINNIillUIMIIMIIIIIIIIIIIklllKlllllllllllllliyilWllUVW 
IOO00O990K0K0-090O0  —  9—  0K9-K0900O0-OK99 
>*.KH*fM*nK9o0K0K99K0r'‘OK09K09-99009O-a9- 
IO—  *^OO0K9N00A.or»o  —  00000-09  KO0-90-K0ON*. 
>N4QAA9><3«MI#AAA49«NA4AAAA0|9O"NAn«All4ANN 


OK  O  —  OO9O0KOO 
9OK90OO9K—  OO 
o-#co--"0»«m 


II  M  M  4  M  W  «U  I 
-O0999O0< 
-  U  9  9  K  0009  i 

pi----- J9| 

•  *•  ••••*•* 
00000904 
»  •  I  I  I  I  I 

9  0  0-00  9  0- 

|  -0  —  00001 

1  -  O  OK  0  o  O  0  4 
B^O-  —  00004 
■KOOOOOOOl 
l  n  *•»•**  • 

H,—  00000«>O4 


>  K  K  0  O  0  - 

>  0  0  O  K  0  i 

IK  0  0  0  0  4 

>9  0  00  0. 

I  0  0  0  M  0  I 

>••«•• 
1000004 

I  I  II  » 

I  —  0  K  O  0  < 

>  9  O  K  K  K  4 

•  090901 
I  (|  O  0  -  O  . 

»  O  O  K  •  •  I 

>000004 


*00001 

>oooo< 

I  III  II  II  V  I 
I  —  —  O  o  I 
I  0  0  9  K  4 
»  O  K  O  O  « 
to-  —  —  I 


I  I  I  I  I  I  I  I  I 

100000000 

>00900000 

I  II  II  II  II  II  II  II  II 
>90000—  K  — 
■K0-0KKOO 
•  0-  K0K0K0 


ooooooooooooooc 
i  I  i  T  i  I  i  i  I  i  T  I  I  I  I 


‘  9  O  —  —  < 
I  •  -  O  K  < 
>0-0—1 
>0090 
>  0  O  OK  I 


IOOOKO90O 
‘—0009000 
‘0009000- 
<00090000 
►  9O--O00O 


0  0  0  0  1 
O  O  O  o  I 

II  II  II  U  I 

O  0  0  O  I 

0  0-04 
O  o  o  K  I 
0  0  O  0  i 

•  •  •  • 

0  0  0  0  4 
I  I  I  I 

O  0  K  0  I 
0  0  0-1 
K  0  •  O  4 


0  0-9 
9  9  0  9 
O0  O  O 
0  0  —  0 
—  —  —  — 
•  •  •  • 
0  9  0  0 
till 

0  o  o  0 
•  00- 
9  0  —  K 


I  I  I  I  I  I  I  I  I  I  I  I 

009000000000 

0OOKO0O0OO—  9 
OKO900OOOOKO 

—  K  —  —  0  —  OOOOO0 

—  9OKO0O-0--- 

000900000000 
I  I  »  »  I  I  I  »  I  I  I 


009990099900 
OOKOaN»9OK09O 
O  —  —  000O00OK0 

009900  »*»  90  JO  *> 


OVOCO0  COO  —  OOOOCOOUOOC  0  0  0  0  4 


>900000000  CO0C,0O9OC9<iO9oOOf,»9 


NNHVIIIIM 
O  9  9  O  9  —  0 
&OO-0OOO 


S%&!S!!!!!!!!Sl!!,l,!^u,k,u,uw^,u,l>v*:^lll,|,,tfVMVUUluiiiikiiiwiiiiinuiiiiiii!ii«iiiiiii 
99O00O00OOK0OO«^O90O0  0O0O0OOK0OO90O9OO90CKKOOO00  —  90 

O«^^*^9999000000000OOOOOKK^0OOOO99OOO--000OO0C90O9O-0O0KOC- 
—  —  —  —  —  —  —  ——  —  —  —  —  —  ——  —  —  —  —  —  —  —  ——  —  —  —  —  —  —  —  —  —  M0000000000000000000000  O  O 

OOOOOOOOOOOOOOOOOOOOOOOO09OOOOOOOOOOOOOOOOOOO9OOOUO 


<00091 

[!3333 

|  o  o  0  9  K 
t  «  «  O  O  0 

I  0  0  0  0  0 


m k  is?; 

8  "  o  9  o  > 

Is  =ss: 

ft  -j  0  o  O  I 

**•*••• 

900001 

i  i  »  T 


>00001 
10  0  0  0  1 

•  4»  •  •  • 

1  ^  W  1 

lilt 


>09  0  9  I 

■09—04 

•  0  0  O  O  I 

>000*4 
►  •  •  •  • 

>00004 
I  I  I  I 

>  0  K  K  0  I 
I  0  0  O  0  . 
’  0  O  K  0  4 

>  •  •  K  K  I 
>99991 

>  •  •  •  • 

>00004 

>00004 

>  O  K  O  —  4 

>  0  9  O  9  > 
•09  004 
>00—04 


i  **5  = 

>  O  9  K  O 
1—0*0 
*90  0  0 
*0000 

>  *>0  0  .*  < 
I  I  I  I  I 


I  V  N  II  II  I 

>  •  •  0  O  I 

•  0  O  K  9  I 

>  0  O  O  9  « 

•  —  0  0  K  « 

>00001 
i  •  •  •  • 

>  o  C  O  o  I 

I  I  I  I  I 

>0009# 

•  0  —  O  9  • 

I  9  K  0  9  » 

10  0  0  0  4 

>  9  9  9  9  * 

\  •  •  •  • 

>  OK  O  0  i 

>009  Of 


'sssif* 

•  90000 
I  9  —  0  «>  K 


II  II  II  U 
9  ®  O  0 
0  0-9 
9  —  0  • 
0  9  0  0 
0  0  0  0 


0  0  0  9  0  4 
-OO0O« 

0*09—4 
O  O  O  O  K  4 

000004 

I  I  I  I  I 

009004 
—  O  K  0  O  4 
0  0  0  9*4 
0  0  0  0  0  4 
9  9  9  9  9  4 
•  •  •  •  • 
0  0  0  0  0  4 

0  0  —  —  o  4 

o  K  o  -  O  < 

0  0  0  0  0  4 
K  O  O  K  0  4 
0900-4 
——0004 

0  0  0  0  0  4 


II  U  II  w 
9  0  0  0 
•  —  0  — 
•  009 
0  O  O  O 
K  K  K  O 


0  0  O  0 
9  9  9? 


I  II  U  II  II  I 

*  O  9  O  C  : 
»  K  0  9  0  4 

>009  —  4 

»  O  «  «  K  4 
>00004 

r.«  •»  O  W  « 

I  I  I  I  I 

*00001 

>  o  o  e>  ~j  < 

IIIIIUIMlI 

r  o  •  o  0  ( 

«  •  0  9  9  < 
*00004 
>000  04 
1  0  0  0  9  4 

*  O  O  O  O  4 
111*1 

>90  —  0  • 
■  K  0  K  9  4 
*90001 
>99  0  01 


HI  IU  *1  W 
0  0  K 
0  0  0  O 
0  0  0  9 
K  *■  K  A. 

0  0  0  0 
•  •  »  • 

—  *.*•»  .*  < 
(III 


9  O  K  O 
0-00 
9  0  0  0 
9  0  0  9 


O  0  K  0 
0  0  •  0 

9  0  0  0 


o  o  o  —  « 

0  0  0  0  4 
0  0  K  O  * 
9  9  9-< 

»  •  *  i 

0  0  0  0  4 
I  I  I  I 

-009* 
0  9  0  0  1 
K  O  O  •  • 
K  K  «  0  I 
«  O  O  O  1 
•  •  •  • 
0  0  0  0  4 

-0994 
O  K  0  •  I 
K  —  O  -  I 
-  •  O  —  * 
K  K  O  9  * 


•  O  *.■  *•  O  v‘  V  O  i> 

IIIUIIIIUUUII 

>*.00  —  0KI»  — 
>0000  —  0^  0 
•VOWN*  —  00  4» 

*000  0  0000 

•  <•  ■-*••»  0  O  **  '» 

I  I  I  I  I  I  I  I  I 

>90000000 

>  o  j  C"?  4>  «  «  r* 

jumiiiiiuviii 

>400000900 

»  —  09  44I0K9  — 


Figure  8  -  Data  Output  for  Computer  Program 


Figure  8  (Continued) 
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C  PROGRAMMER  N.  L  ILL  I  STUN  •  COOE  74S.  EXT  3357 

C  OaCENS 1  TV  IN  SLUG5/CU.FT,H*0EPTH  IN  F  T  »PO* I NI Tl AL  INTERNAL  PRESSURE 
C  IN  LBS/S0.IN.80*IN|TIAL  RADIUS  IN  I NCHE S .C = SPEED  OP  SOUNOIN  MEDIUM  IN 
.C  IN/SEC • P A*ATMOSPHER IC  PRESSURE  IN  LUS/SO . I N , STF *ST ANDOFF  IN  INCHES 
C  ANO  TIME  IS  IN  MILLISECONDS 
C  INITIAL  VALUES 

NEACIS* 101 )  M  M 
101  FORMAT! 12) 

01 MENS ION  TSTPC  3, 1C0C I.STFPC!  3.  1000) . S TFUC ( 3. 1 000 ) .STF (3) 

DIMENSION  A<  100  >).U(  l'JOO)  «0U(  10 Or  ) 

0  IMENS  ION  HC(33)  .UODISO  l.STFOHSO)  ,STF02(50>  .STFD13.S0)  .POD  (SO) 
READ ( 5 • 99 )(HD(J). SOD! J).STFDl(J). STF 02 (J).(STFO(I .J). 1*1.3), 

IPOO! J). 3*1. MM) 

99  FORMAT (8F10.4) 

001  J* 1 • MM  «  1 
H*HO( J) 

80*800! J) 

STF 1*STF01 ( J ) 

STF2*STF02( J) 

OG  9  1*1.3 
9  STP ( I)*STFO( 1. J) 

PC-POO! J) 

C  VALUES  OF  CONSTANTS  FOR  THE  LIQUID  -  WATER 
0*2.0 
e*4.41E4 
C*8.QE4 
P  A*  14.7 

PL*C*h*  32 . 2/1 44 .0  4P A 
C  FEACING 

wniTE(6.88)  80  .  H.PO  ,STF1  .STF  2",  (  S  TF  (  I  )  ,  I  *  1  .  3  ) 

88  FORMAT! 1HI/1H3/////42X.44H  GILMORES  SECOND  ORDER  APPROXIMATION  FOR 

1  THE///45X, 14H IMPLOSION  OF  A.F4.1.19H  INCH  RADIUS  SPHERE///45X , 1 4H 

2  AT  A  DEPTH  0F.F8.1.14H  FEET  OF  »A TER/ / /33X , 5 1 HWHEN  THE  INITIAL  IN 
3TERNAL  PRESSURE  IN  THE  SPHERE  IS.F5.1.5H  PS I A///47X . 33H  ANO  THE  ST 
4 ANOOFFS  ARE.  IN  INCHES///48X, 1CHSTAN00FF  I .F2 I . 2///48X . 1 OHST ANDOFF 
5  2.F21 • 2///48X . lOHSTANOOFF  3.F21.2///48X.1 OHSTANOOFF  4.F21.2/// 

84 AX . 1C HST ANDOFF  S.F21.2) 

C  CALCULATE  INITIAL  VALUES  OF  THE  RADIUS  AND  VELOCITY 
82*12. G736£4)*( PO-PL ) / ( B0*0*2 .0  I 
83*4. C  4B24  4  2/! 3.C*C >- ( 2 .784eE 4 ) *P0»B2/ ( 80*0 *C ) 

84* 3. 9*B2*83/C-4.C* (82* *2  )/(  3. 0 *00)4(2 . 3736E4)  *  ( -P0*83/ ( BO»D«C ) 
l-PL*82/( 3 . 3  *0*80 • • 2)44, 3»PU* (02**2)/(3. 3*80*0*C**2 I ) 
03*2.a»(B2**3)/(C4BO)-2.9*B2*H3/BO4t.8*!B3**2)/C43.2*B2*n4/C 
l4(2.C73t>)*('2.'3»PL*e3/(O*eO**2)-8.0*P0*B4/(BO*C*O) 
2424.'J«PU«,B3*B2/(B0*D*C**2)  ) 

B6*l  ',8.:  »!92  ••2)«83/!  1 5 .0 *C *80  )  4  4 .  S#B4 *83/C4  l  »  Q4B2*B5/ (  3.04C ) 
1-11.2*824*3/(9. ; 480**2 ) -2 . 9*B 3** 2/ ( 2 .C *80 > -9. 4*82*84/ ( 3. 0*80) 

24! 2.9736F4)  *( 1 . 2*PO*B3» *2/( bO*D »C**2 ) ♦ 6. 4«P0*B2«-B4/ ( 3 . 0*80*D*C**2 ) 

3- 2.0  4P0»B5/( 3.3*80*0*C )-PL*B3/( 1 :  .0*D*BO**3) 

4—  J.2»PL*82**2/( C*BO*«  3)-4.9*PL*B4/( SC. 0*D4B0*«2) ) 
fl*0  •  D 

S* ( 1.44C— 4) • SORT (D)*P04*( 1 . C /3. C ) *BO/PL** ( 5. C/8. 0 > 

•HITEI6. 100 ) 

C  HCAOINGS 

180  FORMAT! 123H1  TIME  RAOIuS  VELOCITY  ACCELERATION  BOVP 


HU  1 


fcFN  SOURCE  STATEMENT 


IFNI S) 


I  TSTFl  STFPCI  STFUCI  TSTP2  STFPC2  STFU 

2C2) 

0C66  1*1.4,  I 

A(  I  )  =  eO  +  fci2*H * *24e3«R» *3*H**2*b4«R**24R * *2.i,85#R*#3  +  R*#3*B6*R* »3 
U{  I  )*2.  .  ■H2»R+3.5*B3»R*«244.0*R*B4*R»*2+'g.C*R4*2»e5»R**24 
I6»O*H**3»06*R*#2 
eCYP*P0*B0««4/A<  I  I •• 4 

OU<  I  )*<  <U(  I  )»»3»*<  (  <PL4B)/(B0YP4B>  )••(  3.C/7.  .*  )  )/( 2.0*0 

1- 3.  :*U(  I  )4»2/2#0-7.C«(PU+B)/(  6.  l*D ) »'(2 .C736E4  )  #(  1.0 

2- <  (UOYP*fl)/|PL40)  MMf  .0/7.  ;)4UI  I  )/C*(  (PL4B )/(BDVP+B>  )*»<3.<V7.0> 

3- U { I >/C«t (BCYP40)/<PL4B>)*»(3.C/7.0) ) 

4- 4.; *u(  I  )»B0YP/(0*C )*l 2»  i7  36E4> •( < { BDYP4H ) / < PL+B > >**(4. 0/7. O’ 

5- U( 1 >/C*<PL-*6  J/IBOYP+B) >)✓( A<  I >  —  A  « l)*U( I J/C*< (PL48)/(B0YP 
64E)  >•*< J.C/7.C  )  ) 

Y*A<  t )»UC I  )• *2/2 .C ♦( 2.073664 1 4A ( I >*(80 YP-PL ) /O* (1.3  — 

1 <2.J736E4)*(eOYP-PL )✓<  2.-*0*C»*2> ) 

YKsU ( I  Ml A(  I  >  **2  )»( C»*3)*< 1 .O-Ull >*«2/(2.C»C*#2) >/Y**2 
1-A|  1  )*(C*»2>*( 1.C-U1 1 >/C)/Y 
C4STFU=Y«C/STF 14 VK*Y»*2/(C**3*STFl»*2) * <C*»2- Y/STF l 
1* ( YK**2/C*»3  »*( Y*44/C«*3)/(2.C»STF1»*4  )  ) 

STFUC1=C4STFU/C««2 

STFPCI =0*1  Y/STFl-STf-UCl«*2/2.U)/<2.0736E4) 

14C*( (Y/STf 1 -STFUCI »«2/2. 0)«*2 )/ (2.0* ( C * »2 ) *  I  2. 0736E4 ) ) 

TSTF1=( ISTFl-AI I)  I/O*  I  l.C-UI I  MAI I)/(C»STF1) ) 4R 
C4STFU=Y»C/STF2+YK*Y»*2/<C**3*STF2»*2>*(C4*2-Y/STF2 
14(YK*«2/C«*3)*C Y**4/C**3)/(2.0*STF2**4 ) ) 

STFUC2=C4STFU/C««2 

STFPC2=0*( Y/STF2— STFUC  2**2/2» 0J/I2.C736E4) 

14C*( (Y/STF2-STFUC2«*2/2.<M*2)/(2.C*CC*»2M(2.0736E4> ) 

TSTF2=I ISTF2-AI  I > ) /C ) *  I  1 . O-UC I ) *AI I)/1C»STF2) MR 
CO  709  KK=1,3 

C4STFU*Y«C/STF(KK )4YK*Y**2/(C*»3«STF{KK)**2)»CC**2-Y/STF(KK> 

14( YK*42/C»»3)*{ Y444/C* *31/(2. v*STF(KK) **4>  ) 

STFUC IKK, I )=C4STFU/C**2 

STFPCtKK, I )=C*C Y/STFIKK l-STFUCIKK, I  M*2/2. 0  )  /  I  2.  G736E4  ) 

MO*  I  (Y/STF(KK)-STFUC(KK,I  )  2/2 .0  M*2)  /( 2. 0* I C**2  M ( 2.C736E4 )  > 

7C9  TSTFIKK, I )*( (STFIKK )-A< 1) )/C)*l l.C-UI I ) *  A | I >/(C*STF(KK>  TMR 
R*R*  1 « E3 
TSTF1=TSTFI*I.0E3 
TSTF2=TSTF2*1.0E3 

•  RITE (6. 1 13 1R.AC I ) ,U( I) ,0U( l ) ,BOYP , TSTF 1 . S TKPC 1 , STFUC 1 , TSTF2  » 
1STFPC2.STFUC2 

1 10  FORMAT!  Fll .6.FB.S.2E13.S.E12.S. 2IF9.6, 2E13.31 > 

H*R«1  .OE-3 
TSTf1*TSTF1»1.0E-3 
TSTF2=TSTF2»1.0E-3 
66  R=R+S 
T*R 

C  INITIAL  VALUES  FOR  PREDICTORS  ANO  CORRECTORS 
UP4=U<4) 

AP4*A{ 4 l 
C4*U (41 
04*  A  1 4  ) 

M*4 

C  USING  THE  FOUR  CALCULATED  initial  values  above,  begin  integration 
00669  MS.KCO.l 


RU  I 


EFN  SOURCE  STATEMENT 


I FN ( S ) 


M*M  +  1 
250  L*0 

700  UP5=U1 i-4 )+4.C*S*t 2.< *0U1 l-l )-DUt I-2)+2.**OU< 1-3) )/3.0 
AP5a  A  (  I-4)+4.0*S*12.0*Ut I- I )-Ot  1-21^2.0‘UI I  — 3  > )/3.C 
230  C00=UP5-1 12.9*1 UP4-C4)/1 21  .0 
00C*AP5-1 12.v*t AP4-C4)/12I  .  : 

DC0D» 1  (  COO* •  3  )  * {  1  1PL  +  B  )/t P0*1  80/000  )**4+B)  )  *  •  (  3.  S /7. 9  )  )  / 1 3,  .C  *C ) 

1- 3. . *COO**2/2.0-7.0*1PL+B )/16.0*O >*12.  C7  36E4)*  <  1 .9 

2- <  ipo*i  eo/oooj  *«4+e )/ipl+b)>**ic.c/7.o  >+cod/c*  < ipl 

3  +  B)/tP0*lB0/D00  >**4  +  B>  ) ** 1 3 . 0 /7 .0 >-C0D /C* 1  <  PO*  ( B0/00D )  **4  +  B ) /  C  PL 

4  +  e) )••!  J. 0/7.0)  ) 

5- 4. : *C00*P0*1 60/000 )»*4/l  0*0 •< 2 . 0736E 4  )  *  1  1 (P0*( BO/OOO) **4 

6  ♦  B  )  / 1 PL 4b) )  **l  4. 0/7. 0  l-COO/C*  1 PL+B)/ IPO* (80/000) **4+B) ) )/<000 
7-COO *C00/C*1 ( PL  +  8 )  /  ( P0# ( BO/OUD >  4*4*0 ) ) • • ( 3. 0/7 . 0 ) ) 

C5*19.3*U1  1-1 )-Ul I-3)+3.J*S*tCCOD+2.0*DUt  1-1 >-0Ul  1-2)  )  )/8.0 
05* 19.0* A (  1-1  )-A(  I-3)  +  3.9*S»<  COO+2. 9*Ut  1-1  )-U(  1-2!  ’  )/8.0 
Q*  AOS ( UP5-C5  ) 

IFtG-l..V)2tC'.2lC»260 
c  half  interval  procedure 

400  S*S/2.C 

A4MH*(8-J  .0*AI  1-1  ) +  135.0* A(  T-2)+4C.C*At  I-3»+A(  1-4)  )/256.0 
4+s*l-15.9*Ul  1-1 ) +90 •  C  *U  1 I-2)+15.0*U( 1-3) >/128.0 
A3MH*l 12.C *A( !— 1)+135. 9*A( 1-2 )♦ 1*8 .9 *A l I- 3 ) ♦ A l 1-4 ) )/256.0 
4+S*l-3.  *U1 1-1 )-S4.i*Ul  I- 2 ) +27 • C *U ( 1-3) ) /I  28.0 
04VHa  (  8*-  .C*Ut  1-1  )  +  135.0*U(  I  -  2  )  +  40 . 0*  U(  I-3)+U(  1-4)  )/2S6.0 
4+S*l-15.G*0Ul 1-1 )+90.C*OUt l-2)+15.0*0U( 1-3) )/128.0 
U3MH*( 1 ?.0*U( I- 1 )  + 135.9 *U( 1-2 )  +  lC8.,*U( 1-3) +U(  1-4) )/256.0 

4  +  S* 1  — 3 •  *OUl 1-1 )-54.C*001 1  —  21+27. 0*OU(  1-3) ) /l 28. f 
A (  1-4  )*A3MH 

A( I— 3)*A( 1-2) 

A ( I— 2)sA4MH 
U(  1-4  )*U3MH 
Ut I-3)*U( 1-2) 

U1 I-2)*U4MH 

BCYPsPOMUO/Al  1-4)  )  •* 4 

0U( 1-4 )=( lut I-4)«*3)*U 1PL+B) /l BDYP+B) )** 13.0/7.0) )/( 2.0*0 

1- 3. v  *U( 1-4 )• *2/2.0- 7.C*t PL+B )/t 6.0*0)* 12.0736E4)*! 1.0 

2- 1 tUCYP  +  ei/tPL  +  f) )**t  6.0/7.3 )+Ut  I-4)/C*l tPL+8)/tBDYP 
3+e) )  *•  1 3.0/7.9)-Ut 1—4 )/C*  1 ( BDYP+B )/( PL+B ) ) **13.0/7,0) ) 

4-4.0 *U(  1-4  )«BDYP/t  0*0*1  2.0  736E4)*! ( 1 DC  VP+B ) / ( PL 

5  +  e) )*»( 4.0/7. C )-Ul 1-4 )/C*lPL+B)/t BDYP+B) )  >/» AC  1-4) 

6—  A 1 I-4)«U1 l-4)/C*ttPL+B>/tBDYP+B) ) **13.0/7.0) ) 

OUt  l-3)sOU(  1-2) 

eCYP*P0«< 80/A< 1-2) )**4 

0U1  1-2 )s< (U(  I-2)**3)*l 1 1 PL  +  B) /I BDYP  +  B)  )** 13. 0/7.0) ) /I  2.0*0 

1 - 3. *  *Ul I— 2)**2/2,0  — 7.C*1 PL+B )  / 1  6  . 0*0 ) • 1 2. 0 736E4 ) *11,0 

2- 1 1 BDYP+B )/irL+B) )**16.0/7.0)+U1 I-2)/C*l 1PL+B)/1BDYP 

3  +  6) )«* 13. C/7.9  >-Ut I - 2 ) /C* l I  00 YP  +B  > / 1  PL  +B ) )**13.0/7.C ) ) 

4-4." *Ut  1-2)  *00YP/ 10*0*1  2.0  736E4)*»  1  1  BDYP  +  B) /l  PL 
5  +  B )  )*•( 4.3/7.: )-Ul  1-2  )/C*t PL  +  B) /t BDYP+B) ) ) /I  At  1-2) 
b-Al  1-2)*U1  I-2)/C*l  1PL+B)/1B0YP  +  B>  )  **  1  3 . 0/7.  >7  )  ) 

L=L  + 1 

IFtL-i: )7C0. 700.71 
C  SET  UP  FOR  ITERATION 

266  UP5  =  C5  +  V • C • 1 UPS-C 5)/)21.0 
AP5*D5  +V» 3*1 AP5-05 ) / 1 2 1 .0 
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IFN(S) 


L>Lil 

IF(L-2>230»40C.4C0 
C  CALCUL*T  ION  OF  FINAL  VALUES 
210  U( I >=CS*9.C*(UP5-C5  >/12J .0 
500  A(  I)  =  D54V.')«(AP5-DS)/121.0 
BCYP»PO*BO**4/A<  I  )  4*4 

0U( I  >»( IU( 1 >**3>*< ( (PL4B)/(B0YP+£1) >»*( 3.C/7.C ) >/(2.C*C> 

1- 3.  :  »U(  (  )»*2/2.i>-7.0*(PL  +  B)/(  6.«.*O)*(2.0?36E4)  #(1.0 

2- (  (8DYP*8)/(PL-»B) )**( 6.0/7. 0)+U( I > /CM  ( PL ♦  «>/( B0YP4B >->•** ( 3.0/7 .0  I 

3- U( 1 )/C»{ (BOYP*B)/(PL+B> )•»( 3. 0/7.0) ) 

4- 4.. »U( I ) *BOYP/( 0*C)*( 2. )736E4)*( ( (BOYP+B )/( PL*B >>**( 4.0/7. 0 ) 

5-  U(  I >/C*(PL+fi>/(BOYP+B>  »>/(A<  1 )-A( 1)*U< I ) /C*( (PL*B)/(BDYP 
6  +  L) )•*( 3. 3/7.0)  ) 

V  =  A( I) »U(  I )*»2/2.0*<  2.C736E4) *A« I ) * ( BO YP-PL ) /O* ( 1 . 0- 
1(2.  ’736E4 )*<BOYP-PL  >/( 2.‘  *0*C»*2> ) 

Yk=U( I )»( A( I ) **2 )*(C**3)*( 1 ,0-U ( I )**2/( 2.0*C**2) >/Y**2 
1  —  A (  I )*(C**2)«<  l.w-U(  I  )/C)/Y 
C4STFU=V«C/STF1*YK*Y**2/(C**3*STFI*»2>  * (C**2-Y/STF1 
!♦ ( YK**2/C**3 ) *( Y**4/C**3 )/( 2.C*STF1**4)  ) 

STFUC1=C4STFU/C**2 

STFPCl=C*( Y/STF1-STFUC1 **2/2.0 >/( 2.C736E4) 

1*C*( (Y/STF1-STFUC 1«* 2/2.0  >**2 )/(2.C*(C**2)* (2.0736E4) > 

TSTFlal (STFl-AI 1 ) >/C)*( 1.0- U( 1) *A( I ) /( C*STF1 > )*T 
C4STFU* Y*C/STF  2* VK* Y**2/( C**3*STF2**2) * (C**2— Y/STF2 
1  ♦  ( YK.**2/C*  *3 )  *(Y*«4/C**3)/(  2. C*  STF  2**4  )  ) 

STFUC2==C4STFU/C**2 

STFPC2=0*(Y/STF2-STFUC2**2/2.C.)  /  (  2.  3736E4) 

1*C*( (V/STF2-STFUC2*«2/2.0)**2)/(2.0*(C**2)*(2.0736E4) ) 

TSTF2=( l S TF 2— A (  I  ) )/C)*(  1 .0-U(  I ) • A ( 1 )/(C*STF2) ) *T 
DO  7)8  KK *1.3 

C4STFU=V  *C/STF(KK)  ♦ YK • Y**2/  (  C**  3*  STF  IKK)**2)*(C**2— YZSTF(KK) 
1*(YK**2/C**3)*(Y#*4/C**3)/(2«G*STF(KK)**4)  ) 

STFUC (KK . I )=C4STFU/C**2 

STFPC(KK. 1 )*0*( Y/STF(KK >-STFUC(KK.  I ) ** 2/2 . 0 ) / ( 2 . 0736E4 ) 

1+0* ( ( Y/STF (KK) -STFUC (KK, I ) * *2 /2 .0 ) *• 2 ) / ( 2 . 0* ( C ** 2 ) * ( 2 . 0 736E4 ) ) 

708  TSTF ( KK •  I  >*< ( STF ( KK  )-A(  I ) )/C) *(  1.0-UtI  ) * A ( I >/(C*STF(KK) ) )*T 
T=T*1.0E3 
TSTF1=TSTF1*1.0E3 
TSTf2=TSTF2*1.3E3 

»R1TE(6, 169)T.A( l ).U< l ),DU( I ) .BOYP , T STF 1 , STFPC 1 . STFUC 1 . TSTF2 , 
1STFPC2. STFUC2 

169  F0RMAT(F1 1 . 6.F8 .5. 2E 1 3. S.El 2. S* 2( F9.6. 2E 1 3.5 ) ) 

T*T»1 .Ot-3 
TSTF1=TSTF1*1»QE— 3 
T ST F 2= TSTF 2* 1 .BE— 3 

C  RELOCATE  CERTAIN  QUANTITIES  FOR  THE  NEXT  STEP  OF  THE  INTEGRATION 
C4  =  C5 
04  =  05 
UP4=UPS 
AP4aAPS 

I F ( A( I )— A( I- 1) >669,669,240 
240  IF(U( I) )7J»669»669 
669  T=T*S 

GO  TO  7' 

71  MRITEI6, 75 ) 

7S  FORMAT! 57HTHC  PROCESS  IS  NOT  CONVERGING  QUICKLY  ENOUGH  AT  SOME  STE 
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-  EFN  SOURCE  STATEMENT 


tFNIS) 


IP) 

TO  tF!M-nC?)7.223.223 

223  MRITe(6,22A) 

224  FORMAT ( 9* H  1000  VALUES  HAVE  BEEN  CALCULATED  BUT  THESE  VALUES  00  NO 
IT  COMPLETE  THE  ENTIRE  FIRST  PERIOD) 

7  MR t  T€ ( 6  « 137) 

107  FORMAT! 120H1  f STF3  STFPC3  STFUC3  TSTF4 

1  STFPC4  STFUC4  TSTF5  STFPC3  STFUC5  ) 

00  7 1 S  1*1*3 
00  715  M* I « 1000 
715  TSTF!  I.N)«*fSTF!  I.H)*1.0E3 

MfllTE 16.731 ) ( 1TSTF! I,N).STFPC( I *N) , STFUC ( I *N> . 1*1 .3) .N-1.1000) 

701  FORMAT! 3!FIC .6.2E15.6) > 

1  'CONTINUE 
222  STOP 
END 
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